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A simple and efficient algorithm of the molecular-dynamics simulation of the hard disk 
system based on the Event-Driven method is developed. From the analysis of algorithm, 
the complexity is O(logTV) per 1 event, and the constant coefficient of the complexity 
is smaller than conventional efficient algorithm based on the concept of Cell-Crossing 
Event. The maximum performance of more than 460 millions of collisions per CPU-hour 
on the Alpha600 compatible in a 2500 particle system is achieved. An extension to the 
infinite-space system based on this algorithm is also proposed. 
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1. Introduction 

The molecular dynamics method (MD) was worked out in a paper for the first 
time by Alder and Wainwright in 1957,0 and they subsequently wrote a series of 
works.cHlfl Their simulations were performed with many hard disks (or hard spheres 
in 3D) near the liquid-solid phase-transition point, and they found that the system 
was crystallized despite the particle had only repulsive force. These discoveries over- 
turned the common opinion of those days, and greatly influenced the development 
of the study in computer simulations. 

In the hard disk system, the dynamics consists of only collisions and straight-line 
movement. Since distinct events occur one after another in time, we do not need 
to integrate the differential equations with a constant time step based on Newton's 
equation of motion. The method that is based on the finite constant time step and 
integration with the equations of particles step by step in time is sometimes called 
"Time-Step-Driven Molecular Dynamics" (TSDMD). On the other hand, in the 
hard disk system, the simulation that proceeds based on events is called "Event- 
Driven Molecular Dynamics" (EDMD). Compared with TSDMD simulation, the 
algorithm of EDMD simulation is completely different. We need the knowledge of 
an advanced algorithm and a data structure to perform the efficient simulation in 
EDMD. The strategy of direct computation of particle-pairs result in the complexity 
0{N 2 ) for large particle number N . The point of improvement in the speed in hard 
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disk system is how well we deal with the queue of a future event and the data 
structure. 

The improvement of complexity in the algorithm of large-scale hard disk system 
was developed by Rapaport.!'! This algorithm is based on the concept of sub-cell 
method]] and both Collision Event and Cell-Crossing Event are stored into Multiple^ 
Event Time List. Then the minimum time is searched by Binary Search Tree(BST)J3 
When the event occurs, the particle-pair or the particle — sub-cell respectively 
relevant to Collision Event or Cell-Crossing Event is deleted, and collision time 
for the particle relevant to the event is re-computed and new nodes are created. 
The BST is reconstructed by inserting these new nodes, and the minimum time is 
searched. On this algorithm, the averaged complexity per event becomes O(logiV), 
and the reduction of a large amount of computation is realized. 

However, since the algorithm of Rapaport is very complicated and difficult to un- 
derstand, several algorithms to simplify a data structure and improve the effici ency 
in the large-scale molecular dynamics simulation were proposed in the 90s.lH13I1j 
Marin et al.E3 developed the technique of Local Minima Algorithm (LMA) to avoid 
additional re-computation for Event List. When we actually schedule future event 
list, LMA put only the minimum event time relevant to each particle into Com- 
plete Binary Tree (CBT). In 1995, Marin and Corderdi3 compared various C(log N) 
searching algorithms actually in EDMD simulation, systematically. They concluded 
that the efficiency of Complete Binary Tree (CBT) was the most suitable for hard 
disks system in all density regions. If the number of particles increases, for CBT it 
was clearly shown that efficiency increased significantly when compared with other 
searching algorithms. If one compare CBT with BST, their complexity are the same, 
i.e., O(logiV); but the algorithm of CBT is much simpler than that of BST, conse- 
quently CBT requires less memory space and realizes better actual performance. 

In this paper, we developed an algorithm based on a strategy different from that 
of Cell-Crossing type. The algorithmis extended to Exclusive Particle Grid Method 
(EPGM) developed by Form et al.Jla (Sec. 2). Then, a bookkeeping methodta is 
applied (Sec. 3). Compared with thfi-Cell- Crossing type, our algorithm extended 
the concept of Linked-Cell MethodtZHiS and Neighbor List, which are often used 
in TSDMD to carry out an efficient simulationJlj From the analysis of complexity, 
we show our algorithm is smaller than the complexity of Cell-Crossing type. By 
an empirical evaluation of the simulation in hard disk systems, our code could be 
shown to perform better than that of any previously-published works. 

The infinite volume extension of the sub-cell method is usually considered impos- 
sible because the required memory is proportional to the volume in the conventional 
sub-cell method.L3 We developed the method of compressing information about the 
infinite sub-cells into limited finite arrays. In addition, the hashing method, which 
is considered the most efficient searching algorithm, is applied to our method in 
order to pull out the information on a neighbor cell in high speed. It is found that 
the hashing method can be applied especially easily on our method. 

Various applications in a very wide-ranging field will be possible by changing 
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the external field and the collision rule in the large-scale EDMD. Typical examples 
performed by EDMD so far are: phase transitioninthc equilibrium system (solid- 
liquid transition and two-dimensional melting), IU'ErHo the nonequilibrium fluid sys- 
tem (e.j* Rayleigh-Benard convection), E30 the nonequilibrium chemjcal-reaction 
system JH3 the n o n cquilibrium dissipative system (granular system) P 3 H 24 H 2 ^ random 
disk packingrjO When large-scale computation becomes possible in these systems, 
the value of simulation with a hard disk system must increase significantly. 

The outline of this paper is as follows: in Sees. 2 and 3, the algorithms devel- 
oped are explained. The comparisons with the algorithm of Cell-Crossing type by 
analyzing the complexities are shown in Sec. 4. The empirical evaluation is given 
in Sec. 5. The extension to the infinite system based on developed algorithms is 
explained in Sec. 6. In Sec. 7, a short summary and a final comment are presented. 

2. Extended Exclusive Particle Grid Method 

The sub-cell method is often used to achieve the significant reduction of computa- 
tion. In the EDMD, the complexity of producing and updating future event list for 
restoring event time of each particle-pair are reduced to O(N) and 0(1), respec- 
tively, by the sub-cell method. LCM is the method of dividing the system into small 
sub-cells. When the size of the sub-cell is bigger than a particle diameter, we will 
have a difficulty in coding a program, because we do not know how many particles 
enter in each sub-cell. Therefore, link lists must be prepared for the particles in 
each sub-cell. 

On the other hand, another efficient sub-cell method, called Exclusive Particle 
Grid Method (EPGM), was independently developed by Buchholtz and PoschelLL-H 
and Form, Ito, and Kohring 113 to simulate the soft-core granular particle system in 
TSDMD. In this method, there is only one particle in each sub-cell. Though the 
EPGM is essentially the method of putting a particle to one sub-cell in LCM, it 
does not need to use pointers for neighbor sub-cells or link list. Here sub-cells in 
EPGM are called "grid" . In the EPGM, the length of grid l gx is determined by the 
following inequality: 



O < lgx < V2a, (1) 

where a is the radius of particle. One example of the number of grids n gx , n gy in 
the system of length l x ,l y and the length of grid l gx ,l gy can be calculated by the 
following equations. 



l gx = INT(l x /(V2a)) + 1 
l gy = INT(l y /(V2a)) + 1 



(2) 
(3) 
(4) 
(5) 
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Figure 1: A typical example of the mapping pattern using EPGM when packing 
fraction v = 0.70. Both the position of hard disks and the occupied lattices are 
shown. 



Note that the total number of grids is n g — n gx x n gy . As an analogy to the lattice 
spin system, EPGM is regarded as N + 1-states Potts model of the square lattice 
system. This is because the particles are completely mapped into each lattice (i.e., 
one grid corresponds to one particle respectively), and we put into the rest of 
grids in which there is no particle (Fig. . Since continuous and random positions 
of the particles are mapped into the lattice, the specification of neighbor particles 
becomes very easy. 

Form et al.t-3 applied this algorithm in the high-density soft-core granular system 
with short-distance interaction in TSDMD. They achieved high efficiency on a vector 
computer. Based on this algorithm, the extension of EPGM to EDMD in the hard 
disk system is developed. 

When the system is in high-density, EPGM can be simply applied to EDMD. For 
a candidate of the next colliding particle-pairs, we need to search only 24 neighbor 
grids, which form the square mask. If neighbor grid is not 0, the collision times of 
candidates of colliding particle-pairs are computed only by the registered particle 
number in the square mask. We call these 24 neighbor grids MIN, because this is the 
minimum mask in the simulation of EDMD. Note that if the smaller mask is used, 
the computation will break down during the simulation, since a possibility of overlap 
between a central particle in the mask and particles outside the mask occurs. When 
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Figure 2: Minimum mask (MIN) and larger masks, which are generated by the 
algorithm of making circle on the computer display, are shown. The solid lines in 
each mask denote the minimum distance from the central lattice of mask to the 
frame of mask. 



EPGM is applied in the high-density system, the computation is optimized because 
a sufficient number of particles contained in the mask MIN. These are enough to be 
candidates of particle-pairs of collision, and only the required particles are registered 
in the mask MIN; the efficiency increases as a result of the computational reduction 
of collision time for neighbor particle-pairs. 

On the other hand, when the system is of low density, no sufficient number of 
particles as candidates of collision are registered in the mask MIN. Under such a 
situation, the computation will break down, since collision time between the central 
particle in the mask and particles out of the mask will become the next minimum 
time. In order to prevent this breakdown, the extension of EPGM is developed. 
The region must be extended to look for candidates of collision particle-pairs bigger 
than MIN. Since the rigorous isotropy of neighbors in EDMD is not necessarily 
demanded, the shape of the mask might be approximated by a rough circle. It is 
found that the mask approximated by the lattice-like grid with the circle can use the 
algorithm describing the circle on the discrete space of computer display. Figure ^| 
displayed the circles from R= 3 to R= 10 on the discrete space of computer display, 
which MIN is also showed over Fig. 0. The total number of neighbor grids (mask) 
are 24 (MIN), 36 (R=3), 68 (R=4), 96 (R=5), 136 (R=6), 176 (R=7), 224 (R=8), 
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292 (R=9), 348 (R=10), respectively. This extension of EPGM is called Extended 
Exclusive Particle Grid Method (EEPGM). 

Compared with LCM, EEPGM is simple, because rough neighbor particles can 
be simply regulated by grids so that only the necessary minimum may be taken. In 
EEPGM, since link list and pointers for neighbor sub-cell is not necessary, the re- 
quired memory is also sharply reduced, the number of operations for setting EEPGM 
is small, and the program becomes very simple. Moreover, the extension to the in- 
finite volume system is easy when using a hashing method to EEPGM. This is 
explained in detail in Sec. 6. 

3. Neighbor List and Dynamical Upper Time Cut-Off 

Though EEPGM produces a significant reduction of computation compared with 
considering all particle-pairs, it is inadequate when large number of particle simu- 
lation are performed. In this section, the next step in the strategy of the increase 
in efficiency based on EEPGM is developed. Here, the concept of Neighbor List 
(NL)E3 is adopted. Since the grids correspond to each particle, we can regard the 
registration of Neighbor List as already being completed. Therefore, we need only 
to search neighbor particles along in the form of a mask. 

In the usual way of LCM+NL, after the system is divided into sub-cells, particles 
are listed into the link list. Then neighbor particles within radius t^l are entered 
into Neighbor List from the link list. However, since the length of both lists is 
unknown, they must be estimated by trial and error. Although registration of 
NL is completed by one operation in EEPGM, LCM+NL must use two different 
unknown size lists, which is accompanied by a complicated procedure and requires 
difficult programming. Therefore, EEPGM (+NL) is simple, which means that 
both debugging and extension do not require excessive effort; moreover, only the 
minimum nearest particles can be seen, because the system is divided into the 
minimum sub-cells, i.e., the grid. Since registration of NL is completed by one 
operation in EEPGM, efficiency is better than LCM+NL at a result. 

The next improvement in speed is that the computation of collision time only 
from particle-pairs which are registered in the mask of EEPGM during the time 
t^L- Then the complexity of EEPGM for every event is reducible to 0(1) instead 
of O(N). After time ijvio the grids are again re-mapped in order to update neighbor 
particles. The time of Neighbor List ijvx must be determined such that the central 
particle does not collide with a particle completely out of the mask. If i/vz, is too long 
and such a event happens, the program picks up a wrong pairs of particles to collide, 
which could result in negative collision time. This conventional determination of 
fjvx needs a lot of trial and error. In order to overcome these difficulties, t^L 
is determined by the following procedure. First, after completing EEPGM, the 
maximum velocity v max in the system is searched, and the value of its velocity is 
restored (the complexity of this searching is the same order of EEPGM O(N)). 
Next, time tjsrLmin of the system is calculated. In this calculation we suppose both 
the central particle and the particle out of the mask have the maximum velocity and 
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Table 1: The minimum distance for each mask. 



those particles undergo head-on collision. If t NL — t NLmin , a counting mistake of 
collision pairs in the system never occurred during the time tNLmin- The minimum 
NL distance rjvimm is required when tNLmin is calculated, which becomes clear 
from the geometry of the adopted mask shown in Fig. |[ Therefore, tNLmin is given 

by 

, _ TNLmin ^^max 

tNLmin = r (v) 

^max 

where <J ma x is the maximum radius of particle in the system. The minimum dis- 
tances r m i n for each mask are shown in Table [l]. 

When we simulate in the equilibrium system, this strategy will work well because 
tNLmin hardly changes. However, in the nonequilibrium system (e.g., the relaxing 
process or the system with heat bath) it breaks down because the maximum velocity 
changes drastically at every step. To overcome this difficulty, we must check the 
maximum velocity for each event with energy increase. Fortunately the complexity 
of this checking process is 0(1). Therefore, in the simulation of the nonequilibrium 
system t^L will change one after another. We call this change in t^L techniques 
Dynamical Upper Time Cut-off (DUTC). The development of DUTC, EEPGM 
became applicable in the nonequilibrium system. 

Although we do not need to update the grid pattern for a long time in the 
high-density system, we should often do so when we use the mask MIN in low- 
density systems because the grid pattern changes drastically. In order to reduce the 
frequency of updating grid pattern, we can only use a bigger mask. 

4. Analysis of Complexity 

Analysis of complexity is one of important factor in estimating the efficiency of 
algorithms. In this section, a comparison of analysis of complexity between the 
algorithm of the Cell-Crossing type and the EEPGM + DUTC is shown. The 
difference between the algorithm of Cell-Crossing type and the strategy of EEPGM 
+ DUTC is Cell-Crossing Event itself. Therefore, especially with regard this point, 
both complexities with a constant coefficient k are estimated in the case of A x N 
collisions being actually simulated. Note that the particle number N is supposed 
to be quite a large number and the techniques of improvement in the speed are also 
used in both algorithms described by Marin et al.,E3 

• Cell-Crossing type (LCM + Cell-Crossing Event) 



— The initial and last step (xl) 
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Linked-Cell Method — k LC M x N 

Computation of Event Time — kpp x 9 x N c x N + k PC x4xiV 
Construction of Complete Binary Tree — kcBT x N 
Update the final position of particles — kupp,ATE x N 

- Iteration Step (loop) 

( A x N) x (kpp x 9 x N c +k C BT x log N) + (B x N) x (k PC x 'i+k C BT x log N) 

where N c is the number of particles per sub-cell. The most important point 
is that the Cell-Crossing Event occurs at a certain ratio to Collision Events. 
Therefore, the additional computation of B x N times Cell-Crossing Events 
arc needed when we want to simulate A x N times Collision Events. Since 
the complexity of the terms related to Cell-Crossing Events is always O(N), 
this is not negligible in the actual simulation. 

• EEPGM + DUTC 

- Update of EEPGM (xC) 
EEPGM — k E EPGM x N 

Computation of Event Time — kpp x N g x N 
Search the Maximum Velocity — kuvs x N 
Construction of Complete Binary Tree — kcBT x N 
Update the final position of particles — kjjpDATE x N, 

- Iteration Step (loop) 

(Ax N) x (kpp x N g + k C BT x log TV) 

where N g is the averaged particle number of the mask. Actually the value of 
C is an order C ~ A, which is the same order of the frequency of updating 
EEPGM. 

Now, the comparison with order N in both algorithm with constant coefficient is 
as follows: 

• Cell-Crossing type 

(kLCM + kpp x 9 x N c x (A + 1) + k PC x (3 x B + 1) + k C BT + k UPDATE ) x 
N + (Ax k C BT + B x k C BT) x N x log N 

• EEPGM + DUTC 

(k E EPGM + kpp x2xN g + k M vs + k C BT + kuPDATE) xCxN+(Ax k C BT) x 
N x log N 

The most striking point is that the complexity of Cell-Crossing Events is of order 
of 0(N log N) (i.e., B x kcBT x N x log N). This result of analysis suggest that the 
efficiency of EEPGM + DUTC is better than for Cell-Crossing when the simulation 
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with an enormous number of particle is performed. On the other hand, in the 
comparatively small particle system, the coefficient of C x N terms in EEPGM 
+ DUTC may be larger than Cell-Crossing type. However, the difference might 
be quite small, and it is impossible to estimate an exact coefficient of algorithms 
analytically. To the author's knowledge, the coefficient of N terms is strongly 
dependent on the ratio of N c to N g , and the rough estimation shows both algorithms 
are same when N g /N c ~ 4. Actually, the increase in computing the Cell-Crossing 
Events has almost no effect to the CPU time as the simulation is actually performed 
in a system with a very large number particles. 

5. Empirical Evaluation 

When actually running the simulation, the order of the complexity is not so reliable, 
because it depends strongly on a constant coefficient when the number of the particle 
is relatively small. Moreover, the efficiency of the code changes significantly by the 
ability of a computer, the language, the performance of compiler, and the ability of 
programmers. Though a perfect comparison of efficiency of codes developed by the 
past workers is impossible, some publications showed how many particle collisions 
per CPU-hour can their codes be computed by their computers. 

Marin et al.El simulated with hard disk system, and they also compared their 
code with the codes based on two main high-speed algorithms proposed by Rapaport 
i and Lubachevsky,i respectively. As a result, it was shown that the efficiency 
of the code of Marin et al. was higher for the entire density region. Therefore, 
the code proposed here should be compared with the code of Marin et al., and 
is computed with equal number of particles (N — 2,500). Marin et al. achieved 
the maximum performance of 16.07 millions of collisions per one CPU hour on a 
SUN690 workstation. Note that only the highest performance is shown because 
it is different for different densities of the systems. On the other hand, the code 
proposed here achieved a maximum performance of 460 million collisions per one 
CPU hour (Alpha600 compatible, DEC-Fortran) . It was found that a high efficiency 
was realized even if the performance of the machine was reduced. 

Note that the workstation of our laboratory could actually simulate a 2,500,000 
particle system. In this simulation, the amount of installed memory was 250 
megabytes and the computation performance was 210 million collisions per one 
CPU hour. 

6. Extension to Infinite Volume System 

In this section, a simple example in the open boundary system that does not have a 
ceiling is considered. This is the case when there is an energy source at the bottom 
of the system under uniform gravity. 

The system is also divided into grids by EEPGM. However, because the top of 
the system is open, the grid goes to the top of the system infinitely. This means that 
the number of arrays for the grid becomes infinity, and the simulation is impossible 
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from finite memory. To overcome this difficulty, the hashing method, which is well 
known as the fastest searching algorithm, is applied to keep the number of arrays 
in finite size and to simulate the dynamics of the system with high-efficiency. 

• Construction of Data Structure 

Firstly, the serial number of grid NG(i g ,j g ) is defined by 

N G = N gx (j g -l)+i g , (7) 

where N gx and i g (\ < i g < N gx ) are the total number of grids and the index 
of grid in the horizontal direction, respectively; j g (l < j g < oo) refers to the 
index of grid in the vertical direction. In addition, the maximum number of 
the serial grid Ncmax is calculated by N Gmax = N gx (j gmax - 1) + i gmax ; the 
pair (i gm ax, j g max) is at the maximum grid pair of containing the particle. 

The serial grid Ng(1 < Nq < Nomax) is a one-dimensional array, in which 
particle number or is listed. This is called the one-dimensional Virtual 
Array. When j max is large value, there are many O's in the one-dimensional 
Virtual Array. Though there is only information that a particle is not just 
in grid, this relates to the memory capacity directly. Now Virtual Array is 
compressed. After all, the only information necessary is the particle number 
and its index of grid. Therefore, one-dimensional integer arrays are prepared 
for A(N), BX(N), BY(N), and grids of in Virtual Array are ignored, and 
then packed in order from the end; A(N) stores the particle numbers only, and 
BX(N) and BY(N) stores indexes of i g and j g for each particle, respectively. 
If you want to know whether there is a particle in grid (i' g ,j g ), you need only 
search the index-pairs correspond to (i' g ,j g ) in the lists of BX(N) 7 BY{N) 
linearly However, this process is inefficient because the complexity of 0{N) 
is necessary in searching. 

Therefore, the hashing method, otherwise known as an algorithm which can 
realize the searching with 0(1) is applied. The following simplest hashing 
function is explained here as an example though various hashing functions 
can be considered and there is still room for development. 

First, a hashing function is defined by 

Nr - 1 

k = INT{— ) + f (8) 

which means that the Virtual Array is equally divided by the length L (e.g., 
5 ~ 10) and key k is calculated corresponding to serial number Nq. Then, 
kmax — INT(NGmax/L) is calculated using Ncmax- Then additional arrays 
C{k max ) 1 D{k max ) are prepared, these arrays are restored in C(k) where k 
begins in A(N) and in D(k) how big the arrays for each k. In this case, 
necessary arrays are A(N), BX(N), BY(N), C(k max ) 7 D(k max ), and these are 
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confined to a finite value. Since additional arrays needed to use the hashing 
method are only C(k max ), D{k max ) (k max < N), the amount of memory used 
is only slightly increased in comparison with the linear searching. 

• Searching Process 

In order to know what the particle number is in the grid (i'g,j' g ), the following 
process is carried out. First, a serial number is calculated by N{j = N gx (j' g — 
1) + i' g . Next, k! is calculated by the equation of hashing function (||). Then, 
the searching ranges of BX, BY are limited only to index of [C(k') ~ C(k') + 
D(k') — 1] (< L). If there are equal pairs of grids (i' g ,j' g ) as a result of 
the searching of BX,BY, the index s of BX,BY reveals A(s), which is the 
particle number in the grid (i g ,j' g )- When equal pair is not found in BX, BY, 
there is no particle in the grid. 

This procedure becomes possible in high-speed simulation. The complexity 
becomes 0(1) instead of the linear search O(N), because a search is only 
carried out on the length of L with hashing method. 

A computation is possible for other boundary conditions with the same proce- 
dure if one-dimensional Virtual Array can be created. Therefore, high-speed simu- 
lation is possible in principle for any kind of boundary to be applied. There is room 
for improvement in hashing functions, which is the easiest when dividing equally, 
because it is obviously inappropriate when particles are distributed heterogeneously. 

The strategy of EEPGM has an advantage that it is easily extended. The ease 
of a development is an important factor of the evaluation of an algorithm. 

One problem in low-density systems is the overwhelming increase of the arrays 
assigned to grid in comparison to the number of the particles. This is not desirable 
as the memory capacity is the same of as that of infinite system. However, the arrays 
for the grid are compressed to the size of the particle number in the same way as 
described above. Since supplementary arrays are made by the hashing method and 
information on neighboring grids is efficiently obtained, there is no problem for both 
efficiency and memory capacity. 

7. Concluding remarks 

In this paper, we developed an algorithm for a hard disk system without using 
Cell-Crossing Event, which is simple, efficient and easy to extend. EEPGM can be 
easy to extend because of its simplicity, which can never be realized in LCM. One 
example is the system that hard disks with various size of diameters coexist. Though 
there was a limitation inthe degree of the poly-dispersion with EPGM described to 
Buchholtz and Poschel, LL.3 EEPGM can be applied easily to those systems. First, 
the grid is made based on the smallest particle radius in the system. Next, we have 
only to check the nearest grids by using a suitable mask of the bigger level when the 
poly-dispersity increases. This way, EEPGM has a wider application than EPGM. 
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This code achieved the fastest simulation speed in the world; about 460 million 
of collisions per CPU hour for the 2500 disk system on the VT-Alpha-600. Since 
the order of complexity is O(logiV), the increase of complexity is slow when the 
particle number increases. Now, we can carry out large-scale molecular dynamics 
simulation (<~ 10 6 ) on the usual Workstation in our laboratory. 

Hard particle simulation is a powerful tool for studying the fluidized state de- 
scribed by the kinetic theory or hydrodynamics. Therefore, the large-scale simula- 
tion in the hard particle system will become an important subject. 

Finally, the algorithm in this paper is suitable for the scalar machine, and the 
development of an algorithm for the parallel machine is the subject of a future 
study. Note that we discuss only hard disk systems in 2D for simplicity, but an 
extension to hard spheres system in 3D is easy to be carried out. 
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